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Perturbative calculations in field theory at finite temperature involve sums over the Matsubara 
frequencies. Besides the usual difficulties that appear in perturbative computations, these sums 
give rise to some new obstacles that are carefully analized here. I present a fast and realible 
recipe to work out sums over the Matsubara frequencies. As this algorithm leads to deal with 
very cumbersome algebraic expressions, it has been written for computers by using the symbolic 
manipulation program Mathematica. It is also shown this algorithm to be self-consistent when 
it is applied to more than one loop computations. 



*This work was partially supported by the CICyT (Spain) under contract AEN/93/673. 
f Address after January 1, 1994: Department of Physics and Astronomy, Northwestern University, Evanston, 
Illinois 60208, USA. 



A system of particles at temperature T = 1/(3 is completly described by its partition 
function. From it, it is very easy to get the other physical magnitudes, e.g. pressure, energy, 
entropy, etc. So, one of the main tasks of a statistical field theory is to develop methods to 
obtain the partition function. This can be acomplished by means of a perturbative expansion in 
the coupling constant. Following this idea, one finds that the partition function is represented 
by the one-particle- irreducible Feynman diagrams without external legs (bubble diagrams). 
This is just a simple example to point out the importance of perturbative methods in finite 
temperature field theory, which have been widely and sucessfully applied to different purposes 
as the QCD deconfinement. 

At finite temperature, together with the usual difficulties that a perturbative computation 
involves (integration over loops, regularization and renormalization), some others arise because 
of the sums over the Matsubara frequencies. The time-like component of the momenta is 
quantized because the definition of the partition function requires the fields to be periodic in 
time (fermions are anti-periodic). Its allowed values are the so-called Matsubara frequencies: 
i2mr/j3 for bosons and i(2n + l)n/(3 for fermions (n G Z). This makes the integration over 
that component change into a sum over frequencies 

We will focus on evaluating sums over the Matsubara frequencies. The method that I will 
describe in what follows is well known but I think is worth noticing its self-consistency on 
evaluating each sum over the Matsubara frequencies of more-than-one loop diagrams. Also, a 
computer program which provides a fast and realible way of performing these sums is presented. 
Finally, some examples are worked out by using such program. 

Let us sum a function / over the set of frequencies k° = iu n = i 2nn / (3 



f J nGZ 

If / does not have any singularity along the imaginary axis, we will multiply it by a function 
with simple poles and residue one at k° = iu n , e.g., l/(2iti)(/3/2) coth(fc°/3/2). So, (|l|) is written 



where Up{z) = |coth(/?z/2) and C is the contour in fig. |T](a). This contour can be deformed 
into the one shown in fig. |l](b) and, therefore, eq. @ can be expressed as 




(1) 




(2) 





Using the contours in fig. [TJ(c), we can write the integrals in (|3|) as 




and 
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Now, if f{z)up{z) goes fast enough to zero as \z\ — > oo, the integrals at infinity will vanish. In 
this case, we conclude 

Wf{k G = i 2nn/(3) = - £ Res {f(z)u p (z)} , (4) 

where z a are the poles of f(z). 

A similar result can be obtained for fermions at zero chemical potential 

^T,f( k ° = *(2n + 1)*//?) = - £ Res {f(z)v p (z)} , (5) 

where vp{z) = \ tanh(/3z/2). 

Considering the dependence on the momenta in a Feynman diagram we can go further in 
our analysis. The contribution of an L-loop Feynman diagram with 4-momentum kj = (k®, kj) 
flowing around each loop (j — 1, . . . , L) can be written as 

where = i uj n] is the corresponding Matsubara frequency, / also depends on k± . . . k^ and in 
general has the form 

/(*?, . . . , k\ ) = £ h (k°i, ■ ■ ■ , hi) n h <J«h° + a R , b R ) ; (7) 

R=l V i J 

The D-function is defined as 

£>(a;,y) = — L- - (8) 
x 2 — y z 

it is a non-analytical function that arises by noticing that, with E\ = k 2 + m 2 , D{k G ^Ek) is 
the propagator for a scalar field, (i/k + m)D(k° , Ej,) is the one for fermions, —g^ u D(k°, \k\) + 
(1 — ^)knk p D 2 (k°, \k\) is the one for photons and similarly for other fields. The function h is 
analytic on k®, . . . , k° L ; of = {+1, 0, —1}; a R and b R do not depend on and the exponent 7^ 
is an integer. The global sum means that / may be a sum of functions with these properties. 
An interesting property of the Feynman diagrams is that if a R , of, of and af 7^ 

o?of = of of VR * R>. (9) 

In other words, if a propagator in the diagram carries momentum k — p + E 1; there may exist 
other propagators with momentum —k + p + E 2 , but not k + p + E 3 nor —A; + p + E 4 , i.e. the 
relative sign of two momenta, e.g. k and p, is constant. 

A detailed and straightforward analysis shows that the result of performing a sum over, 
let us say, the jth Matsubara frequency of a function that verifies the previous properties is 



another function that verifies the same properties. So, if we were able to find a method to 
evaluate such a sum, we would be able to perform iteratively every sum. 
Let us use eq. to compute 

i £/(*?,...,*£) (io) 

for a scalar field. Using the ([?D and @ one finds that when af; 7^ each /^-function has two 
poles of order 7^ at 

fc£ = -^fE^*? + «?±& H ) =4- (ii) 

It is worth stressing that the function / has no other poles but the ones given by the .D-functions. 
So, 

i £/(*?,...,*£) = - E (r«{/«^)}+ Res {/«,(*)}) 

= E Res + (6* - -ft*). (12) 



R 



Where, 



1 rllR- 1 r 1 



The result for fermions would have been obtained simply by changing up{z) into vp{z). Equa- 
tion (|13|) and its fermionic equivalent are all we need to perform a sum over a single frequency; 
additional sums can similarly be computed by following the previous analysis. 

So far, we have not considered whether evaluating these sums is completly sensible without 
mentioning the possible divergences. We just evaluated them by using a regularization method; 
so they do not indeed diverge. The expected divergences that will arise when calculating the 
space-like integrals should be handled by the usual renormalization schemes. 

Summarizing, the recipe consists on identifying the poles of the .D-functions and evaluating 
the residues at these poles. Then every residue is summed. The result is another function on 
the remaining momenta. 

However, this apparently simple scheme cannot be easily acomplished for diagrams with 
a large number of propagators and loops. The amount of algebra nedeed when dealing with 
physically interesting cases become an akward task to be performed by hand. Computers, 
though, have shown to be good at doing cumbersome and repetitive tasks. 

Symbolic manipulation programs have been developed to handle algebra. Some of them 
are intended to perform perturbative calculations in physics, e.g. Schoonship 0, REDUCE [[J 
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and FORM M. On the other hand, there are general purpose programs, like MAPLE [||] or 
Mathematica |J, that provide functions that allow one to do additional manipulations on the 
results, e.g. integrate over the space-like momenta, plot results, etc. In this context, I present 
here a Mathematica package that evaluates the sums over the Matsubara frequencies and that, 
therefore, helps to calculate Feynman diagrams at finite temperature. I have called this package 
FreSum^. 

Let us compute the one-loop correction to the self-enegy of a scalar field interacting with 
itself through a A</> 4 potential. It will be (see fig. |2|) 

A r d 3 k 1 v 1 

" 2 J (2vr)3 (3 % (kP)* - El 

= ^r^e^Dik^) (i4) 

Once Mathematica is started, the following commands are needed to perform the sum over the 
quantized time-like momentum: 

In[l]= «fresum' 

In [2]= Even[kO] ; 

In [3]= f[kO_] = dd[kO,Ek,l] ; 

In [4]= resultO = FreSum [f [kO] ,k0] 

The first line instructs Mathematica to load the package FreSum. The following instruction 
defines k° = i2nir/f3 as an even Matsubara frequency (if we had considered a spinor field, it 
would have been an odd Matsubara frequency: i{2n + l)vr//5, and we would have used the 
FreSum function Odd[]). The third line defines the function f(k°) that we want to sum; here 
dd[] is the Mathematica .D-function. Finally we call resultO to the result of performing the 
sum over k°; it is 

up{E k ) 
Ek 

Now, we can go further and take the massless limit 

In [5]= result 1 = resultO /. Ek->k; 

We can put aside the dependence on the temperature with the FreSum functions VacuumPart [] 
and MatterPart [] ; the vacuum part is the T — ► limit of the argument, while the matter part 
is the result of taking away the vacuum part from the argument. 

lr This package is available at the directory pub/math/paquetes of the Anonymous FTP site ftp.uam.es. 
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In [6]= result2 = MatterPart [result 1] 

In [7]= result3 = result2 /. Nb[x_,bb] -> l/(Exp[bb x]-l) 

we get —Ng(k)/k, where bb is (3 = 1/T, k = \k\ and 

^M S ^*H = S55£FT (15) 

is the bosonic occupation number. If spinor fields are involved, the result will be written by 



using 



1 1 



mz)s .- vM = — mTl (16) 

which is the fermionic occupation number. The integration over the 3-momentum can be 
performed: 

In [8]= result4 = lambda k~2 result3 / (4 Pi~2) 
In[9]= result = Integrate [result4,{k, , Infinity}] 



Finally, the result for the self-energy is 



which coincides which the one obtained in |T| and for this case. 

Another interesting calculation that can easily be done with the help of FreSum is the 1- 
loop photon self-energy with zero external momentum. The 00-component can be written in 
the Feynman gauge (£ = 1) 



IT° (0) = -4e 2 J d 3 k - J2(k 2 + El)D\k , E k ). 

" fc 



The instructions 



In[l]= «fresum' 
In [2]= Even[pO] ; 
In [3]= Odd[kO] ; 

In[4]= f[kO_] = (kO"2+Ek"2) dd[kO,Ek,2] ; 

In [5]= resultO = -4 ee2 FreSum [f [kO] ,k0] 

In[6]= resultl = 4 Pi k"2 resultO / (2 Pi) "3 

In [7]= result2 = Simplify [MatterPart [resultl] ] 
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give rise to 

n°°(0) = — jf dkk 2 N^(E k ) [N*(E k ) - l] . (19) 

Note that the third argument of dd[] is the power of the corresponding D-function. The 
sixth instruction simulates the angular integration. In this case, the vacuum part vanishes 
in the zero external momentum limit due to the Lorentz invariance of the theory; so, here, 
the MatterPart [] function just changes Uo[] into Nf [] . Now, writing Np(z) explicitly and 
performing the changes of variable k(3 — > k we find that the photon self-energy is 

n» = ^jM (20) 

where the function 



6 [Qo k 2 expVk 2 + ~ 

9( z ) = -5 dk ~ f x2 ( 21 ) 

71 Jo ( 1 + exp Vk 2 + z 2 



is shown in fig. |3|. The following chain of instructions will produce such result: 

In [8]= result3 = Simplify [result2 /. Nf [x_ ,bb] ->1/ (Exp [bb x]+l)] 
In [9]= result4 [k_ ,z_] = Simplify [result3 

/. {bb->l, ee2->l,Ek->Sqrt[k~2+z~2]}] 
In [10]= g[z_] := 3 Integrate [result4 [k,z] , {k, 0, Infinity}] 
In [11]= Plot[g[z] , {z.0,10}] 

As g(0) = 1, we recover the well known result in the massless limit |], || 



n?(o) = ^- (22) 

So far, perturbative calculations have almost exclusively been done by using numerical 
methods. However, when one is concerned with topics that require a non-limit behavior, addi- 
tional algebra arises and, therefore, symbolic manipulation programs are required in order to 
use computers efficiently. This is the context in which the package that we have described here 
will be a useful tool. 

We have seen that the evaluation of the sums over the Matsubara frequencies can easily be 
accomplished. Two simple examples have been worked out to show how the package works. A 
NeXTstation, spent less than 2 sec. to get (|17D . The photon self-energy (|T9"D calculation took 
2.5 sec, while fig. |3] took 1,500 sec. because the required integrations have to be worked hard 
by Mathematica. We can learn a number of facts from these computations; the sums over the 
Matsubara frequencies are algebraically awkward, but they are performed systematically, no 
matter the number of loops. One eventually gets a result that has to be integrated over the 



3-momenta; this turns out to be the main obstacle to completly evaluate an involved diagram. 
This task, when dealing with more-than-one loop diagrams, will probably require numerical 
methods of integration to be achived. 
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List of FreSum Objects. 



Fresum[f , kO] Performs the sum of f over the kO Matsubara frequency. 

£ E fe o f(k°) 

VacuumPart [expr] The T = limit of expr. 

MatterPart [expr] The result of taking away VacuumPart [expr] from expr. 

dd [x , y , N] The Mathematica version of the Nth power of the 

D-function: D(x,y) = l/(x 2 - y 2 ). 

bb The Mathematica version of (5 — 1/T, the inverse temperature. 

Ue[z,bb] The Mathematica version of up(z) = ±coth(/3z/2). 

Uo[z,bb] The Mathematica version of V p(z) = |tanh(/3z/2). 

Nb[z,bb] The bosonic occupation number: Np{z) = l/(exp(/3z) — 1). 

Nf [z,bb] The fermionic occupation number: Np(z) = l/(exp((3z) + 1). 

Even[xl,x2, . . .] xl,x2, . . . are defined as even Matsubara frequencies. 

x n = iu) n = i 2nn/ f3 

0dd[xl,x2, . . .] xl,x2, . . . are defined as odd Matsubara frequencies. 

x n = iu n = i(2n + l)ir//3 
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